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1.   INTRODUCTION 

Visual  texture  is  known  to  play  an  important  role  in  the  field  of 
pattern  recognition.   There  are  three  problems  which  are  under  study  in  the 
context  of  automated  image  processing,  viz.  texture  discrimination/recogni- 
tion, texture  analysis  and  texture  synthesis.   In  the  literature  many  methods 
are  proposed  to  deal  with  the  above  problems  [7,  10,  15,  IT,  21,  25,  26,  27, 
281.   Most  of  these,  however,  have  little  or  nothing  in  common  and  lack 

generality. 

In  this  paper  we  intend  to  propose  a  general  method  for  the  synthesis 
of  texture.   This  method  is  based  on  a  mathematical  model  which  treats  pixels 
(picture  elements)  of  a  digitized  textural  scene  as  a  two-way  seasonal  time- 
series  . 

1.1    Texture  Synthesis:   Its  importance 

The  synthesis  of  natural  looking  textures  deserves  considerably  more 
attention  than  it  has  received.   General  methods  are  needed  by  which  we  can 
generate  a  scene  that  bears  an  acceptable  resemblance  to  the  texture  to  be 
imitated  ('parent'  texture).   It  is  highly  desirable  that  the  parameters 
required  for  the  synthesis  are  derived  from  the  analysis  of  the  'parent' 
texture. 

Texture  synthesis  procedures  of  this  type  can  be  evolved  to  solve  the 
ever-pressing  problem  of  image  storage.   For,  it  requires  considerably  less 
storage  to  save  a  few  parameters  and  boundary  conditions  than  to  store  a 
complete  digitized  scene.   This  effectively  means  redundancy  in  the  image 
has  been  reduced,  leading  to  a  more  efficient  analysis  of  the  resultant  scene. 
Using  the  values  of  the  texture  parameters  as  features  for  discrimination/ 
recognition  is  also  suggested.   By  specifying  intervals  on  values  of  each 
parameter,  we  can  potentially  generate  a  family  of  textures. 


1. 2  Texture  Synthesis:   Earlier  work 

We  find  that  many  authors  proposed  methods  [7,  8,  13,  20,  26]  that 
are  "based  on  the  subpattern /placement  rule  model  for  texture.   Among  them 
only  Rosenfeld,  et.al  [26]  and  Conroy  [7]  attempt  to  synthesize  natural- 
looking  textures.   While  others  merely  generate  scenes  with  pre-specified 
statistical  properties  to  be  used  in  their  experiments  concerning  the  visual 
perception  of  texture.   In  their  papers,  Rosenfeld  and  Conroy  demonstrate 
that  a  rich  variety  of  textures  can  be  generated  using  only  simple  subpat- 
terns  and  placement  rules.   In  these  methods  the  success  of  synthesizing  a 
natural-looking  texture  depends  strongly  on  the  choice  of  the  subpatterns  and 
the  appropriate  placement  rules.   These  important  choices  normally  are  not 
based  on  the  anlaysis  of  the  'parent'  texture,  instead  involve  from  human 
judgment  on  a  trial  and  error  "basis.   This  is  the  main  drawback  of  these 
schemes  based  on  the  subpattern/placement  rule  model.   It  is  difficult  if 
not  an  impossible  task  even  to  extract  the  subpatterns  (whose  size  and 
shape  are  unknown)  from  the  analysis  of  any  textural  scene,  not  to  speak  of 
decoding  the  rules  of  their  placement. 

1. 3  Seasonal  Time  Series  Model 

The  subpattern/placement  rule  model,  though  easy  to  visualize,  creates 
many  practical  difficulties.   Here  we  adopt  a  different  approach:   Our 
scheme  for  texture  synthesis  is  based  on  a  model  that  views  the  pixels 
of  a  digitized  two-dimensional  textural  scene  as  a  two-way  time  series.   We 
will  then  resort  to  the  well-developed  methods  of  the  time  series  analysis. 
By  analyzing  the  'parent'  texture  we  identify  the  process  that  generates 
the  time  series,  and  estimate  its  parameters.   Then  with  proper  initial 
conditions  we  can  generate  the  time  series  which  in  essence  resemble  the 
pixels  of  the  synthesized  textural  scene. 


In  the  subsequent  chapters  the  time  series  strategy  is  discussed  in 
detail.   The  effectiveness  of  this  method  is  illustrated  by  examples. 


2.   SEASONAL  TIME  SERIES  MODEL  FOR  TEXTURE 

Any  digitized  visual  scene  can  be  viewed  as  a  two-way  time  series. 
The  orderly  repetitiveness  of  a  subpattern  which  is  an  essential  charac- 
teristic of  texture  hints  at  the  possible  representation  of  a  digitized 
textural  scene  by  Seasonal  Time  Series .   The  choice  seems  quite  natural 
and  very  appropriate  because  there  is  a  striking  correspondence  between 
some  of  the  problems  that  occur  in  textures  and  those  that  are  considered 
in  seasonal  time  series  analysis.   For  example,  in  TV  scan  of  a  texture 
the  repeating  subpattern  is  not  necessarily  identical  from  line  to  line 
though  it  retains  similar  characteristics.   In  the  time  series  analysis 
the  minor  variations  of  the  function  from  one  season  to  another  is 
accounted  for,  by  the  assumption  of  the  presence  of  white  noise.   Well 

developed  methods  are  available  to  estimate  the  values  of  the  parameters 

2 
{•u,  a  }  of  this  noise. 

It  is  possible  that  there  can  be  more  than  one  periodicity  present 
in  the  textural  scene,  such  as,  the  subpatterns  consisting  of  sub-subpatterns 
and  so  on.   This  may  be  treated  as  multiple  seasonality.   In  the  case  of 
statistical  textures  these  latter  periodicities  may  be  absent  and  they  may  be 
treated  as  an  ordinary  time  series  with  no  seasonal  effect. 

Spectral  analysis,  in  the  frequency-domain,  comprises  one  class  of 
techniques  for  time  series  analysis.   Our  interest  here  is  to  stay  in  the 
time-domain  and  build  stochastic  models  for  time  series.   This  way  we  can 
gain  a  better  insight  into  the  nature  of  the  system  that  generates  the  time 
series  and  the  models  can  be  used  to  obtain  forecasts  of  the  future  values 
of  the  time  series.   Our  objective,  from  the  engineering  view  point  is  to 
obtain  models  that  possess  maximum  simplicity  and  the  minimum  no.  of  parameters 


consonant  with  representational  adequacy.   Precisely  the  same  approach  is 
taken  by  Box  and  Jenkins  while  dealing  with  this  subject  in  their  book 
"Time  Series  Analysis  -  Forecasting  and  Control"  (Holden-Day,  1970). 

In  their  book  they  discuss  various  stochastic  models  for  time  series. 
They  also  present  methods  for  identifying,  fitting  and  checking  the  models 
for  a  given  time  series.   They  also  extend  their  methods  to  deal  with  sea- 
sonal time  series.   We  have  presently  implemented  these  methods  in  our 
texture  analysis  programs.   In  order  to  preserve  continuity  we  will  briefly 
review  some  of  the  models  and  methods,  while  omitting  much  of  the  detail 
which  can  be  obtained  from  their  book. 

2.1     Stationary  Models 

In  these  models  the  process  that  generates  the  time  series  is 
assumed  to  be  in  equilibrium  about  a  constant  mean  level. 

Let  us  introduce  some  notation  and  operators  that  are  used  exten- 
sively.  Let 

t-29  t-1'  t'  t+1 

be  the  discrete  time  series  which  can  be  denoted  as  [Z  ]. 

A  series  of  values  of  a   ("shocks")  are  assumed  to  be  generated  from 

2 
a  white  noise  process  with  mean  '0'  and  and  variance  a  . 

a 

'B'  is  the  "backward  shift"  operator  such  that 

B  Zx  =  Z^  .  :   hence   B  mZ4  =  Z^ 
t    t-1  t    t-m 

'  V  is  the  "backward  difference"  operator 

v  \  -  \  -  \-i 
•  (1-B)  zt  • 


2.1.1   Auto  Regressive  Models 

In  this  model,  the  current  value  of  the  process  is  expressed  as  a 
finite,  linear  aggregate  of  previous  values  of  the  process  and  a  'shock'  a 

Zt  =  *1  Zt-l+*2  Zt-2+  •'•  +  *1  Zt-p+  \ 
This  is  called  an  ' aut or egress ive  (AR)  process'  of  order  'p'.   Let 

$p  (B)  =  (l-i^B  -  (J)2B2  ...  -  c(,pBP) 


t 


Then 


!p(b)  zt  "  at 


<j>  (b)  is  known  as  'autoregressive  operator'  of  order  p.   There  are 
p 

2 
p+1  unknown  parameters,  <$>  ,  $    ,  cj>   ...  <j>  ,  a   vhich  m  practice  have  to 

be  estimated  from  the  data. 


2.1.2   Moving  Average  Models 

In  this  model  the  present  value  Z  is  assumed  to  be  linearly  depen- 


dant  on  'q'  of  previous  a's. 

Z  =  a  —  9  a    —  y _  a   _  ...  —  y     a 

t     t     1   t-1     2   t-2        q  t-q 

=  9  (B)  a, 
q-     t 

Q  (B)  is  the  "moving  average  (MA)  operator'  of  order  q. 

4. 


2.1.3   Mixed  Models 


In  many  practical  cases  the  parsimony  in  number  of  operators  is 
achieved  by  mixing  the  autoregressive  and  moving  average  terms  in  the  same 
model. 

VB)  \  -  VB)  at 

2. 2     Nonstationary  Models 

Many  series  actually  encountered  in  practice  are  not  stationary  in 
nature.   But  some  of  them  show  homogeneous  nonstationary  behaviour.   Box 


and  Jenkins  show  that  such  cases  can  "be  represented  by  a  generalized  auto- 
regressive  operator  $(B)  in  which  one  or  more  of  the  zeroes  of  the  poly- 
nomial #(B)  is  unity,  i.e.,  $(B)  =  J  (B)  (l-B)   ,  where  <j)  (B)  is  the 

ir  P 

stationery  autoregressive  operator. 

Thus  a  general  model  which  can  represent  homogeneous  nonstationary 

behaviour,  is  of  the  form 

i  (B)  (1-B)d  Z      =   0  (B)  a+ 
P  t     q     t 


4(B)  wt  =  5(B)  a 


where  w =  V   Z 

Homogeneous  nonstationary  behaviour  can  therefore  be  represented  by 

th 
a  model  which  calls  for  the  d  difference  of  the  process  to  be  stationary. 

The  above  representation  is  a  powerful  model  for  describing  sta- 
tionary and  nonstationary  time  series  and  is  called  an  'autoregressive 
integrated  moving  average  (ARIMA)  process'  of  order  (p,  d,  q).   There  are 
(p  +  q  +  2)  unknown  parameters  to  be  estimated  from  the  data. 
0  (B)  i    _1(B)  Sd 
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Figure  2.1.   Block  Diagram  for  ARIMA  Process 


It  can  be  seen  that  the  ARIMA  process  may  be  generated  from  white 
noise  a  by  means  of  three  filtering  operations  as  shown  in  Figure  2.1. 


2.3 


Seasonal  Model 


It  is  shown L5J  that  the  time  series  which  show  seasonal  effect  can 
be  represented  by  a  general  multiplicative  model: 

~*p(B)  ;p(BS)  Vd  VsD  Zt  =  9q(B)  6Q(BS)  at 


where  s   =  Seasonal  period 

Vs  =  (1-BS) 

tj)  (B  ),  6  (B  )  are  the  polynomials  in  B"  of  the  order  P,  Q  respec- 
tively. 

This  is  a  modified  version  of  ARIMA  process,  said  to  be  of  the 
order  (p,  d,  q)  x  (P,  D,  Q)  .   This  multiplicative  model  is  very  useful  in 
that  it  can  be  easily  extended  to  take  care  of  multiple  seasonalities. 

2.h  Identification  and  Fitting  of  Models 

Identification  methods  are  rough  procedures  applied  to  a  set  of  data 
to  indicate  the  kind  of  representational  model  which  is  worthy  of  further 
investigation.   By  'fitting'  we  mean  the  estimation  of  the  various  parameters. 
The  selection  of  the  right  kind  of  model  results  from  the  extensive  inter- 
action between  theory  and  practice.   The  necessary  theory  is  developed  in 
Chapters  2-5  and  the  methods  for  identifying  and  fitting  models  for  the 
actual  time  series  have  been  explained  and  illustrated  with  examples  in 
Chapters  6-9   in  the  Box  and  Jenkins  book.   J 

The  programs  that  have  implemented  these  methods  have  been  described 
in  detail  in  Appendix  A. 

It  can  be  observed  that  the  extensive  use  is  made  of  the  behaviour 
of  the  autocorrelation  function  (ACF),  <j>  ,  and  the  partial  autocorrelation 
function  (PACF),  <j>   in  the  process  of  identification  of  the  proper  model. 
Tables  6.1  (pp.  If 6)      [   Appendix  A9.1  (pp.  329) ^  and  Table  3.2  (pp.  32) ^ 
are  used  as  reference  tables.   The  initial  values  of  the  parameters  can  also 
be  estimated  with  the  help  of  these  tables. 

Chapter  7  of  the  book'-''-'  describes  a  procedure  to  obtain  a  final  set 
of  parameters  for  a  general  ARIMA  process  from  the  initial  values  by 
minimizing  the  'sum  of  squares  function'  S(<j>,6), 


wh  ere 

n  ,  2 

s(<i>,a)   =     Z      [a  1^,0, w] 

t=l     ^ 

<f>   and  8    are   vectors    of  length  p    and  q  that   represent   the 
coefficient  values    of  the  AR   and  MA  operators,    respectively,    and 

W     =  vd  Z 

2 .  5  Forecasting 

One    of  the  primary   objectives    in  "building  stochastic   models 
for   discrete   time   series    that   occur  in   many  practical   cases    is   to  utilize 
their   ability   to    forecast   the    future   values    of  the   series.      These    fore- 
casts  provide    a  basis    for   a)    economic   and  business   planning,   b)    production 
planning,    c)    inventory   and  production    control,    d)    control   and  optimization 
of  industrial  processes. 

Let   Z    (2.  )  ,  2=1,2,...   be   the    function   that  provides   the   forecasts 

Li 

at   origin    t   for   all    future    lead  times    (l  )  .      This  will  be   known    as    the 

'forecast    function'    at   origin   t.      For  a   general  ARIMA   model    (non-seasonal) 

$    (B)    Vd  Z+    =    6    (B)    a+  2.5.1 

p  t  q  t 

the   forecast   function   Z   (2)    may  be   expressed  directly   in   terms    of  the 

difference   equation  by: 

\^]  •  h+i  -  +i  Wi  +---+  Vd  zt«-P-a  -  ei  at«-i 

-■-   Vt«-q+   V*  S-5-2 

Box  and  Jenkins   show  that  the  minimum  mean   square    error   forecast 
Z,  (O    (2    >   0)    of  Z  is    the   conditional  expectation 

V"  ■  izt«]  =  E[zt«IV  zt-i---] 
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The   conditional  expectations   of  the   terms    in   2.5.2    are   evaluated  by- 
inserting  actual   Z's  when   these    are   known,    forecasted  Z's    for   future 
values,    actual   a's   when  these    are  known,    and   zeroes    for   future   a's 
(because   E|a        1   =   0    for  j    >   0). 

This    forecast   function   is  highly  useful  when   the   future  values 
are   needed  for  very   short   lead  times    as    required  in   many  business    and 
industrial   applications.      For   longer   lead  times   the    forecast   error  will 
be    cumulatively  larger.      Also,    as  the   future  values   of  a's    are   replaced 
by   zeroes,  the   eventual   forecast    function   takes   on   a  deterministic 
mathematical   form  as    dictated  by   the   solution    of   the  homogeneous 
difference  equation    containing  autoregressive   terms    only. 

Our  intention    of  fitting  a  time   series    model   for  the  pixels 
from  a  textural  scene   is   to  be    able   to   forecast   the  series    for  longer 
lead  times    and  thus   effect  the  texture   synthesis    from  a  set   of 
boundary   conditions.      The    forecast    function  will  be  unsuitable    for  this 
purpose    for   the   reasons    mentioned  above.      Especially   the  eventual 
disappearance   of  the  stochastic  effect   is  highly  detrimental  for  the 
textural  property  of  the  scene. 

2.5.1      Generation    of  Time   Series 

As   shown   in   Figure   2.1,   the   general  ARIMA  process    can  be 
generated  from  a  white  noise  process  with   appropriate   filtering 
operations.      Figure  2.2   shows   the  rearrangement  of  the   filters   in 
Figure  2.1   for  the    generation   of  time   series.      a     is   the   output   of 
the  white  noise    generator.      The  present  scheme   requires    (p+q+d)    delays 
and  registers   to   store   the  previous   occurrences   of  Z's    and  a's  which 
are   required  in   th   calculation  of  the  present  value.      The  multiplicative 
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Figure   ..,   .       Generation   of  Time  Series    [Z,  J 
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constants  (<t>'s  and  9's)  are  the  parameters  of  the  model.  To  start  with, 
the  registers  are  loaded  with  the  'initial'  conditions  and  the  values  of 
Z  (l    >    0)    are   successively   regenerated. 

The    above    discussion    can  be  extended  to   seasonal  time   series. 


2.6  Synthesis   of  2-D  Texture 

So   far  we   have    considered  one-dimensional  time   series.      By 
concatenating  either  successive   rows   or  successive   columns,    a  two-way 
series    can  be  treated  as    a  one-way  time   series.      Of  course,  by   doing 
this  we   are   introducing  one  more    (known)   periodicity. 

The  pixels    from  a   digitized  textural   scene   to  be   imitated 
('parent'    texture)    are   row    (or    column)    concatenated  to   form  a  one-way 
time   series    [Z    ].      With  the  help   of  programs   USID*,   LIKESURF*    an 
appropriate   model  is    fitted  and  the  values    of  the   parameters    are 
estimated.      The   generation  process   shown  in   Figure   2.2   is   simulated 
in   GENTEX*  which   synthesizes  the  textural  scene   from  a  set   of  boundary 
conditions.      This   method  is    illustrated  with   an   example   in   the 
next   chapter. 

2. 7  Redundancy  Reduction 

The   digitized  TV  scan   of  a  picture  is    a  row-concatenated  one- 
way  time   series.      When    an    appropriate   model  is    fitted  for   this    series, 
we   come  up  with   an   attractive   scheme   for  redundancy  reduction   in  the 
transmission   of  the  TV  scan   of  the  picture.      In  the   literature    [30]  we 
find  a  method  known  by  the  name    'optimum  linear  predictor'    in  which  the 


*     See  Appendix 
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next   value   in   the   series    is   predicted  by  expressing  it    as    a  linear 
function   of  the  previous   occurrences: 

k 

v    =    £     e.  y    . 

*t        .    _      i  Jt-i 
i=l 

e  =  y  -  y   where  y  =  the  predicted  value 

~G  ~C  T»  T> 

y     =  the   actual  value 
x< 

e     =  the   error  in  prediction 

2 
0's    are   selected  to  minimize   the   mean   square    error    [a    (e    )].      Data 

compression    (redundancy   reduction)    is    achieved  by  transmitting  the 

error  values   only.      At   the    receiving  end  the   actual  values    (y's)    can 

be    reconstructed  from  the    'errors'    (e's)    from  a  system  that   is   to  be 

initialized  to   start  with. 

In  the  present  scheme  that  involves   the   time   series   model, 

a  one-step-ahead  forecast   is   made  by  expressing  the  next  value   as 

a  linear   function   of  not   only  the  previous  values   in  the  series,  but 

also   the  previous    'errors': 


Z^  =   E   d>.  Z^_    .    +   Z   6.  a^  . 
1    i-1  X   *"1   1  =  1     ±      t"1 


(autoregressive  terms)  (moving  average  terms) 


\  •  tzt  -  zt] 

where  Z  =  one-step-ahead  forecast  of  Z 
Xi  "^ 

a  =  the  'residual'  or  'error' 

The  values    of    <j>'s    and   0's    are   selected  to  minimize   the   mean   square   error 

2 
[a    ].      We   can   see  that  this   is    an  improved  version   of  the    'optimum  linear 
a 

predictor'    method  for  achieving  redundancy  reduction   in   the   transmission   of 
a  TV  scan   of  a  picture.      At  the   receiving  end  we    can  perform  a  similar 
'  filtering' operation   as   shown   in   Figure  2.2   to  reconstruct   Z's    from  a's. 
Here,    of   course,    'white  noise    generator'    is   replaced  by  the  transmission    line, 


Ik 


3.   SYNTHESIS  OF  TENURE:   A  CASE  STUDY 
3.1     Selection  of  an  Appropriate  Model 

The  textural  scenes  that  are  used  in  the  following  examples  are 
taken  from  Brodatz    .   A  digitized  scene  of  cheesecloth  texture  (D-105, 
Brodatz    )  with  lb   graylevels  is  considered  as  the  'parent  texture' 
(Figure  3.l)  in  example  1.   For  the  sake  of  analysis,  pixels  from  a 
32  x  l6  window  is  treated  as  a  two-way  input  time  series.   It  is  row- 
concatenated  (s  =  32)  and  is  given  as  an  input  one-way  time  series 
([Z  J,  t  =  1,512)  to  the  program  USID*.   The  output  of  this  program  con- 
sists of  the  estimated  autocorrelations  for  [W  ]  where  W  =  V   V    Z  : 

~G  ~G      J_     S     u 

(a)   for  the  original  series,  Z 

i.e.  d  =  0   D  =  0  [ACF-00] 

("b)   for  the  series  differenced  once  with  respect  to  "basic 

interval  only,  V   Z 

i.e.  d  =  1  D=0  [ACF-10] 

(c)  for  the  series  differenced  once  with  respect  to  seasonal 
interval  only,  V   Z 

i.e.  d  =  0   D=l  [ACF-01] 

(d)  for  the  series  differenced  once  with  respect  to  "both  "basic 

and  seasonal  intervals,  V_,  V   Z^ 

1st 

i.e.    a  =  1     D  =  1  [ACF-ll] 

Figure   3.3   shows   the   plots    of  the    autocorrelation   function    (ACF) 
and  partial   autocorrelation   function    (PACF)    for  all  the   four  cases  men- 
tioned above.      The  partial  autocorrelation  function  is  more  useful   in 


*  Appendix  A 
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the   identification  of  the   non-seasonal  models   and  hence   is   not   further 
discussed. 

The  autocorrelation  values  in  case  (a)  are  large  and  highly 
periodic  as  might  be  expected.   In  case  (b)  it  is  seen  that  by  dif- 
ferencing with  respect  to  the  basic  interval  reducing  the  correlation 
in  general  while  a  heavy  periodic  component  remains.   We  obtain  station- 
arity  by  differencing  with  respect  to  the  seasonal  interval  as  seen  in 
cases  (c)  and  (d)  where  the  correlation  values  diminish  very  rapidly. 
The  values  of  d  and  D  which  produce  stationarity  in  the  present  case 
are  0  and  1  respectively. 

A  prospective  model  is  selected  for  further  analysis  based  on 
the  information  furnished  by  the  autocorrelation  function  by  referring 
to  the  Table  3-2    and  appendix  A9.1    .   In  practice  we  can  pick  up 
many  tentative  models  and  even  probably  overparameterize  them.   The 
program  LIKESURF*  is  used  to  estimate  the  parameters  of  a  given  model 
from  analyzing  the  input  data.   This  program  is  capable  of  checking 
many  models  at  a  time  and  its  output  indicates  the  adequacy  of  any  model 
as  well  as  any  redundancy  in  the  choice  of  number  of  parameters. 

By  the  procedure  mentioned  above,  we  select  the  following  model 

to  represent  the  given  data: 

(1-(|,1  BS)  Wt  =  (l-A1  B)(l-As  BS)  at  (l) 

where  Wx  =  V   Z^ 
t     s  t 

s  =  length  of  season  (length  of  row) 
3.2    Estimation  of  the  Parameters 

We  need  to  estimate  the  values  of  four  parameters  (cf>  ,  A  ,  A  ,  a  ) 

_L   _L   s   a 

by  analyzing  the  input  data. 
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Figure    3.1.      Cheesecloth  Texture 


Figure    3.2.      Synthesized  Texture 


LAGS 

C.c:  16. u  3L\U  48.0  6U.C  '  8J.0  36.0  112.0  128.0 
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Figure   3.3.      Serial  Correlations   in  Example  1 
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Figure    3.U.      Output   of  LIKESURF  in  Exajnple   1 
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Figure  3.5.   Serial  Correlations  of  Residuals  in  Example  1 
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Figure   3.6.      Forecast  Values    (shown  by   x)    in  Example   1 
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\  =  (l+*}  zt-s  "  *  zt-2s  +  at  "  xi  at-: 


-    X      a,         +   X.,    X      a  ,    v 

s      t-s  1      s      t-s-1  (2) 

With  the   knowledge   of  a    ,    a        ,    a     „,    etc.      We   can  make    a  one-step- 

ahead   forecast : 

\+l  =    (1+*'    Ws   "   *   Zt+l-2s   "   Al   at   "   A2   Vl-s 

(3)** 

The  'residual'  a,  ,  =  ZJ  _,  -  Z.  .'  ,  i.e.  the  residuals  are  considered  on 
t+1    t+1    t+1 

the  one  step  ahead  forecast  errors.   With  appropriate  starting  values 
we  can  see  that  the  values  of  a  can  "be  recursively  calculated  for  a 
given  set  of  parameters. 


3.2.1  A  Method  to  Obtain  Starting  Values  for  Recursive  Calculations 

of  Residuals 

A  procedure  to  obtain  the  starting  values  is  described  in  detail 
and  illustrated  with  an  example  in  Chapter  9  of  the  book  by  Box  and 
Jenkins    .   It  is  as  follows:   For  a  given  set  of  parameters  the  series 
may  be  forecast  backwards  starting  somewhere  in  the  middle  of  the  series. 
To  begin  with  the  unknown  a' s  are  assumed  to  be  zero.   This  introduces 
transients  which  will  hopefully  die  down  by  the  time  Z  ,  Z   ,  Z   etc. 
are  estimated,  provided  the  starting  point  is  chosen  sufficiently  far 
along  the  series.   The  back-forecast  continues  till  we  obtain  sufficient 
number  of  terms  to  be  able  to  forward- fore cast  the  first  value  (Z  ): 

\   ■  d+*>  Vs  "  *  Zl-2s 

Here  we   need  the   values   till   Zn    _      and  the   values    of   a^    for   t    <   1   is 

l-2s  t 

assumed  to  be  zero. 


**  While  forecasting  the  unknown  values  of  a  are  replaced  by  their  ex- 
pected value.   In  this  case  E(a    )  =  0. 
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Now  a,  =  Z.  -  Z. 

Ill 

From  now  on  we  can  recursirely  estimate  the  values  of  residuals 

a„ ,  a„  . . .  a  .   We  define  a  sum  of  squares  function  '  S'  as: 
12      n 


The  program  ' LIKESURF'  plots  the  sum  of  squares  function  over  a 
given  parametric  space  and  ohtains  the  values  of  the  parameters  for 
which  it  is  minimal. 

The  output  of  the  program  is  shown  in  Figure  3.^-.   The  estimated 
values  of  the  parameters  are: 

d>  =  -0-15   Xn  =  0-25  X     =  0-5 

Is 

and  min. sum  of  squares  =  526.538 

2 
So  the  variance  of  the  residuals  is  a   =  1*03, 

a 

3. 3    Checking  the  Adequacy  of  the  Fit 

According  to  the  model  the  residuals  a  are  generated  by  a  white 
noise  process.   If  the  fit  were  to  be  adequate,  the  calculated  residuals 
should  be  uncorrelated.   Figure  3-5  shows  the  correlation  values  of  the 
residuals.   It  can  be  seen  that  the  values  are  fairly  small  and  there  is 
no  periodic  component  either.   Thus,  it  offers  no  significant  departures 
from  randomness  among  residuals,  confirming  the  adequacy  of  the  fit. 

As  mentioned  in  the  previous  chapter  we  can  also  check  the 
adequacy  by  comparing  the  actual  values  to  the  forecast  values  from  any 
origin. 

Forecasting  is  done  on  the  same  principle  as  before  (i.e.  using 

different   equation   approach).      The   values   of  a,     for  t   >   t      .     .      is    set 

■^  t  origin 

to  zero. 
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Figure  3.6  shows  the  actual  values  (solid  line)  and  forecast 

values  (crosses)  from  the  origin  t  =  128.   Considering  the  limited 

number  of  parameters  steering  the  model,  the  forecast  function  follows 

actual  values  rather  closely  even  for  t  >>  t   .  . 

origin. 

3-k.         Generation  of  Texture 

We  can  create  a  white  noise  generator  having  a  mean  =  0  and 

2 
the  variance  a   ,  which  generates  the  residuals  a  .   To  start  the  recur- 

3,  T-» 

sive  procedure  for  the  generation  of  future  values  we  need  at  least  2s 
values  of  Z .   These  can  he  considered  as  the  boundary  conditions. 
Figure  3.2  shows  the  texture  synthesized  in  this  manner.   Here  using 
two  rows  of  length  6k    (s  =  6k)    from  the  original  scene  as  the  starting 
values,  the  future  values  of  the  pixels  have  been  generated  by  the 
procedure  mentioned  above. 

The  effectiveness  of  this  scheme  can  be  tested  by  attempting  to 
fill  the  holes  in  the  'parent'  texture.   The  results  are  shown  in 
Figure  3.7-   It  can  be  seen  that  the  edge  effects  are  clearly  visible. 
This  can  be  expected  in  the  case  of  textures  which  have  an  associated 
structure. 

Let  us  consider  a  different  class  in  the  second  example.   The 
texture  of  a  handmade  paper  (Figure  3.8a)  belongs  to  a  general  class  of 
statistical  textures.   It  has  no  structure.   The  time  series  model  fitted 
for  this  texture  is  of  the  form: 

(1-0. 8BS)  Z  =  (l+0.2  5B)(l-0.5BS)  a 

X  Xi 

where  s  =  seasonal  period  (length  of  the  row) 

°l  ■  3-2T 

Figure   3.8b   and  Figure    3.8c  show  the  results   of  hole-filling  operation 
using  the   above   model.      Here   edge  effects    are   hardly  noticeable. 
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3. 5  Conclusion 

It  has   "been   demonstrated  that   the  natural   looking  textures    can  be 
synthesized  using  a  time    series   model.      This  procedure    appears   to  be   very 
effective  while   dealing  with  the   statistical  textures.      Many  biological 
textures  belong  to   this    general    category. 

We   may    add  here   that   the  motivation   for  this   work  has   been 
provided  by  our   desire   to   obtain   a  procedure   by  which  we    can  make    a   just 
choice   of   the   size   and  shape    of  template    for  the   R.O.C.    analysis    of 

r  -i    O  "I 

texture  .      It  was   our  belief  that   the  number  of   auto  regressive    and 

moving  average   terms    and  the   values   of  the  parameters   in    a  time   series 
model   for   a  given   texture   influence   that    choice.      At   the  present   time, 
work  is   being  done  to   establish   a  link  between  the  present  work   and  the 
R.O.C.    analysis   of  texture. 


2U 


(a) 


(b) 

Figure  3.7.   Filling  the  Holes  in  a  Textural  Scene  (Example  l) 
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(c) 
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APPENDIX  A 
Description    of   Computer  Programs 

Acknowledgment 

The   description,    flowcharts    and  the   listing  of  the   programs 
found  in   Bacon's    thesis    [2]    and  in   the   appendix  of  Box  and  Jenkins 
"book    [5]   have  been   very  helpful  in   the   development    of  the    following 
programs.      Several  modifications   have  been    done   to  Bacon's   programs    to 
suit   our  needs.      Efforts   have  been   made   to   retain    the   same  nomenclature 
and  interpretation   of  several  parameters    in    order   to    conform  as    far 
as   possible   to  the   description    of  Bacon's   programs. 

Program  1 

USID    (Univariate   Stochastic   Model   Identification) 

1.1  General  Description 

This   program  is    developed  to    assist   in   the   identification   stage. 

It    calculates    the   serial   correlations    and  their  standard  errors    for   a 

series,    [  Z ,  ]    and  any    of  its    differences    of  the    form    [Vn    V      Z    1 .      The   main 
t  1st 

program  USID   controls    the   input    data   and  performs    the   required   differencinj 
while  the  subroutine  SER  is   used  for  estimating  the  serial  correlations. 

1.2  Discription  of  Input 

A  total   of  m  cards    are    required  by   program  as    input    for  each 
series,   m  depending  on  the  number  of  observations. 
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Card  1 

FORMAT    (8lU) 
NC  :  The  series   identification  number 

NT  :  The  number  of   data  in   the   series 

KPER  :  The   length   of  the   period   (s) 

NLG  :  If  NLG=0,    the    original   data  is    analyzed 

If  NLG^O,    the   mean   of  the   series    is 

subtracted  from  each   term  Z    . 
NM  :  The  highest    degree    of   differencing  by 

basic   intervals    (i.e.d.  ) 
WK  :  The  highest   degree    of  differencing  by 

seasonal  intervals    (i.e.D.) 
NL  :  The   maximum  lag  of   autocorrelations 

to  be    computed 
WLP  :  The  maximum  lag  of  partial  autocorrelations 

to  be   computed 

Cards   2-m 

These    cards    contain  the  original   series    data.      The   present   input 
format  is    32 Zl.      This   may  be    changed  as   per  the    convenience. 

Last    Card 

After  the  last  card  of  the  final  series,  a  card  should  be  placed 
with  0000  in  NT  field  (columns  5-8)  which  causes  the  program  to  stop. 
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1.3 


Computation 


1 


n 


Mean   of  the    series  ,    z  =  —     E      z^ 

n   t=l     t 


If  NLG^O,    z|    =    (z   -z),    otherwise    z'    =    z 


VT    =   V?  VD   zx    for  t=l,2 
t  1st 


n 


Mean  W  =  -     E      W^ 
n   t=l     t 


The   autocovariance    function    for   lag  k(C    )    is: 

K. 


i 


n-k 


C.     =  —     E       (w    -w)(w        -v] 
k       n  t  t+k 

U      J_ 


The  autocorrelation   function   for  lag  k(r    )    is 

k 

rk  ■  <Vco 


The  standard  error   for  the   estimated  autocorrelation   r     is:       (Bartlett's    approx. 

k 


d(r   )    =  -f-  {l+2(r^   +    r2   +...+    r2)}1/2 

k  /n  12  q 


k   >    q 


The   partial   autocorrelation    function    for   log  I  ( <J>      )    is 

for  I    =   1 


=  «< 


where    $£J   =   *  j 


£-1 
r     =     E      3>  r 

1  T_n  «--l,J         4 -J 


£  -1 

1   -      E 
J=l 


for  £    >  1 


K£-1,J      J 


J=l,2  ,.  .  .1-1 


I       r£-l,2,-J 
The    standard  error   for  the  estimated  partial   autocorrelation    formula  is 


3(4>      )   =   l//n      for  k  >   p 

KK. 
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1. h  Description   of  Output 

For  each   series,    for   all  values    of   d  =   0   through  NM  and  D  =   0 
thorugh  NL,    the   autocorrelation   function   and  partial   autocorrelation 
function  has   "been   plotted  for  the   specified  number   of  lags.      The   output 
also    consists    of  the   standard  errors   of  each    function   at   every  lag. 

Program  2 

LIKESURF 

2  . 1  General  Description 

This   program   computes    a   grid  of  residual  sum  of  squares   over 
a  specified  region   of  parameter  values   and  provides    confirmation  test 
results   on  the   residuals    from  the  best    fitting  model.      This    is    a 
modified  version  of  Bacon's  program  and  accommodates    any  model  of  the   form: 

-   A        RS+2   -   X,   B2S    -    A,      nB2s+1 

s+2  2s  2s+l 

This  would  include   many   commonly  occurring  seasonal  models. 
Provision  has   been  made   to   accept   a  constrained  form  of  equation    (l) 
such   as  : 

71  's    d-VW  "    d-\B)(l-A/)    at 

The  program  locates   the   minimum  sum  of  squares    and  identifies 
the   parameter  values  which  produce   that  minimum.      The   serial  correlations 
of  the   residuals    are    computed  and  plotted.      Subroutine  WES   computes   the 
penodogram  plot   and  also    computes   the  value   of  Bartlett's   statistic 
M/C    [2].      These  values    are   used  for   checking  the  adequacy  of  the   fit   of 
the  model.      For  more   details   the  reader  is   referred  to  Bacon's  thesis. 
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The   random  shocks    'a    '    are    generated  by   a  subroutine    GUASS 
when   the   appropriate   parameters    are   supplied.      Using  this    subroutine, 
the   future   values   of  the   series    are    generated  with  the  help   of  the 
fitted  model. 

The   subroutine   SHWPAX,    treating  the  time   series    as   the  pixels 
of  the   generated  scene,    rearranges    and  quantizes    appropriately   and 
points    a  two-dimensional   gray-level  picture   of  the   resultant   scene  with 
the  help   of  the    'PAX'    subroutines. 

2.2  Description   of  Input   Data 

It  is  possible  to   analyze   several   series   sequentially.      It  is 
also   possible   to   check   on   several  models    for   the   same    series. 


Card  1 


NC 

NT 

NPER 

NPH 

NLG 


F0PMAT(5lU,2FU.l) 

The   series    identification  number 

The  number  of  data  in  the  series 

The   length  of  the  period  s 

The   autoregressive    lag  r  in    (l) 

If  NLG  is   nonzero,    the   mean   value   is 

subtracted  from  each  term  of  the  series; 

otherwise,   the  original  series   is    analyzed 

The  value   of  d   (0   or  l) 

The  value  of  D    (0  or   l) 


Cards   2-m 

Contains   the   original  series    data.      The  present   format   is 
(32Z1).      May  be    changed  according  to  the   requirements. 
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Card,  m+1 


ND 

NPR 
NORG 


NXX 


FORMAT    UlU 

The  number  of  series    data  to  "be   used 

to  estimate   the  parameters 

The  number  of   data  to  he    forecast 

The  value    of  the   origin    from  which 

forecasts   are  made.      The  values   of 

z      for  t=l  through  NORG   are   used   as 

boundary   conditions 

Is   a  code  variable 

NxxfO    for  models   of  general   form   (l) 

Nxx=l   for  models    of  the   form 

V*  VD    (l-A  Br)z   =(l-A   B)(l-A   BS)    a+ 
Is  It  1  s  t 

Nxx=2    for  models    of  the    form 


Nxx=3   for  models   of  the   form 


Vl  Vs    (1-*1Br)zt=(1-A1B-A2B2)(l-AsBS)    at 


(2) 


Vd  VD(l-4)   Br)z  =(1-A   B)(l-A   BS-AQ   B2s)    a+       (3) 
lslt  1  s2s  t 


(U) 


Cards   m+2,   m+3,    n&h 

These    cards    contain  the   ranges   of  the  values   of  parameters   to 

be    considered.      The   range   is   specified  for  each  parameter  by  the   lowest 

value,   the  highest   value   and  the  step   change.     Each   card  specifies   the 

ranges   of  three   parameters  with   the   format    (9F6*3).      Card  m+2   specifies 

the   ranges    for  ch    ,    An    and   A„ .      Card  m+3   for   A    ,    A      ,    and   A      _.       Card 
11  2  s'      s+1  s+2 

rn+U   for   A  A^        ,    X_    ,  _ .      The   greatest  number  of  values    of  any  parameter 

2s        2s+l        2s+2  J 
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that   can  be   tested  in   one    fit   is   l6.      If   a  parameter  is  to  be   restricted 
to   one    quantity,   say   0.0,   the   three   entries    should  be    0.0,   0.0,    1.0. 
The   last   four    cards   only  need  be    repreated   for  subsequent 
models   to  be   tested  for  the   same   series. 

Last   Card  for  Each  Series 


A   card  containing   zeros    in  ND   field  notifies    that    another 
series   is  to  be  processed. 

Final    Card 

At  the   end  of   final   series,    a   card   containing   zeros   in  NT 
field   causes   the   program  to   stop. 

2 . 3  Description   of  Output 

The   output    contains   the  sum  of  squares    function   plotted  over 
the   defined  parametric   space.      The   set   of  parameters   which  produce   the 
minimum  sum  of  squares   is   shown.      The  values   of  the   residuals    for  this 
set   is    also  shown.      The  serial  correlations   of  these   residuals    are 
plotted.      The   forecast   function   is   also   plotted  starting  from  N0RG. 

NPR  values    of  the  series    are    generated  using  the  values   of 
'a   '    generated  from  the   subroutine   GAUSS   starting  from  NORG.      These 
generated  values    are  passed  on  to  the   subroutine  SHWPAX  and  its   output 
contains   the   gray-level  picture    of  the  two-dimensional  scene,    as 
described  earlier. 
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Program  3 
GENTEX 

3.1  General  Description 

This  program  generates  the   future  values   of  a  seasonal 
time   series    from  a  given  model  for  a  specified  boundary   condition. 
In   a  limited  sense,    in  the   present   context   it    generates   the  pixels 
of   a  digitized  texture   scene.      The   model  and  the  values   of   its 
parameters    are   embedded  in  the   program  and   can  be    changed  for 
various    cases . 

3.2  Description  of  Input  Data 

Only  the  necessary   details    of  the   series    and  boundary 
conditions    are   fed  in   as   input    data. 

Card  1 

Contains  NC,  NT,  NPER,  NPH ,  NLG,  U,  V  in  the  same  format  as 
that  of  LIKESURF.  All  parameters  have  same  interpretation,  except  NT 
indicates   the  number  of  data  included  as   boundary   conditions  . 

Cards   2-m 

Contains  boundary   conditions.      Input    format   at  present  is 
(32Z1)   which  may  be   changed  as   per   convenience. 

We  will   illustrate  with   an  example   a  method  of  embedding 
any  model  in  the  program: 

Suppose   the  model  to  be   implemented  is   of  the   form: 

V?  VD    (l-cj)1Br)z  =(1-AR)(1-A   BS)    a+  (2) 

Is  It  1  s  t 

(Values   of  d  and  D   are   restricted  to  be   0   and  l).      Using  this   model,   a 

one-step-ahead  forecast   is: 


3U 


z\    =   d*z^    .    +   D*z^        -   d*D*z^ 
t  t-1  t-s  t-s-1 

+    A    (-z.         -   d*z+         .    -   D*Zj_  +    d*D*z^  J 

1      t-r  t-r-1  t-r-s  t-r-s-1 

-   An    a,    .    -   A     a,        +   A.  A     a,  (2 

1     t-1  s      t-s  1   s      t-s-1 

The  values   of  the  parameters    of  this   model   tj).  ,    An  ,    X    ,   a 

1        1        s        a 

are   specified  in  the  beginning  of  the  program.      The   equation  2'    is 
appropriately   coded  and  placed  in  the   loop   as    a  one -step- ahead  forecast 
function,  which   is  the   only  thing  that    changes    for  every  model. 


3. 3       Description  of  Output 

The   output    consists    of  the   generated  values    of  the  time   series 
As    mentioned  before,    these    are  passed  on   to   the   subroutine   SHWPAX  to 
point   a   gray -level  picture    of  the   resultant   scene. 
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